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Abstract 

This paper is devoted to numerical simulation of a charged particle beam submitted to 
a strong oscillating electric field. For that, we consider a two-scale numerical approach as 
follows: we first recall the two-scale model which is obtained by using two-scale convergence 
techniques; then, we numerically solve this limit model by using a backward semi-lagrangian 
method and we propose a new mesh of the phase space which allows us to simplify the solution 
of the Poisson's equation. Finally, we present some numerical results which have been obtained 
by the new method, and we validate its efficiency through long time simulations. 

AMS subjects: 35B27, 76X05, 65M25 (65Y20). 

Keywords: two-scale convergence, semi-lagrangian method, Vlasov-Poisson model. 

1 Introduction 

Recent papers have proved that the two-scale convergence theory developed by Allaire in [2j and 
Nguetseng in [17J can be used successfully in order to develop numerical methods for solving ODEs 
or PDEs with oscillatory singular perturbations: for example, Frenod, Salvarani and Sonnendriicker 
have developed a two-scale PIC method in [13j for simulations of charged particle beams in a pe- 
riodic focusing channel, and Frenod, Mouton and Sonnendriicker have developed a two-scale finite 
volume method in [9j for solving the weakly compressible ID Euler equations. We can also cite 
the work of Ailliot, Frenod and Monbet in [1] about the simulation of tide oscillation for long term 
drift forecast of objects in the ocean. 

On the other hand, many papers devoted to the numerical simulation of Vlasov-type prob- 
lems involve Particle-In-Cell methods (see Birdsall and Langdon [3]) or Eulerian methods like 
semi-lagrangian schemes (see Sonnendriicker et a/.[19j, Grandgirard et a/. [14], Filbet and Sonnen- 
driicker Cheng and Knorr [5j). Since papers like [9j or [13j can be viewed as parts of a 
work programme which goal is the development of two-scale numerical methods for simulations 
of magnetic confinement fusion, it can be interesting to couple two-scale convergence results on 
Vlasov models such as two-scale models obtained in [12j with a semi-lagrangian scheme. 

However, it is preferable in a first time to develop such a method on a more simple problem 
in order to study its behavior especially in the context of non-smooth solutions. The context of 
charged particle beams in a periodic focusing external field described in [13j offers a relatively 
simple framework for answering these questions. 

We recall that such a phenomenon can be successfully represented by the 3D Vlasov-Maxwell 
system. In the same spirit of [13j, we will consider non-relativistic long and thin beams, so we 
can replace the full three-dimensional Vlasov-Maxwell system by its paraxial approximation. To 
obtain this approximation, we do the following assumptions: 

• the beam has already reached its stationary state, 

• the beam is long and thin, 

• the beam is propagating at constant velocity Vz along the longitudinal axis 

• the beam is sufficiently long so we can neglect the longitudinal self-consistent forces. 
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• the external electric field is supposed to be /-periodic in z and independent of the time, 

• the beam is axisymmetric, 

• the initial distribution /o is concentrated in angular momentum. 

Under all these assumptions, the 3D Vlasov-Maxwell system reduces itself to a 2D Vlasov-Poisson 
system where the longitudinal position coordinate z plays the role of time, and even to a ID 
axisymmetric Vlasov-Poisson system of the form 

f%r,Vr,t = 0) = fo{r,Vr) , 

K (r, t) = - - i^o r + i^i (cji ^) r , 

where r > is the radial component of the position vector in the transverse plane to the prop- 
agation direction, G M is the projection of the transverse velocity in the transverse plan to 
the propagation direction, e is the ratio between the characteristic transverse radius of the beam 
and the characteristic longitudinal length of the beam, = f^{r,Vr^t) is the distribution func- 
tion of the particles, = E^{r^t) is the radial part of the transverse self-consistent electric field, 
= S^(r, t) is the radial part of the transverse external electric field defined with the dimension- 
less real constant coi and the tension functions Hq and Hi which are respectively constant and 
periodic. All these quantities and variables are dimensionless. This system is naturally defined 
for r > but we can extend it to r G M by using the conventions f^{—r^—Vr^t) = f^{r^Vr^t), 
E'^{-r,t) = -E',{r,t), and S^(-r,t) = -S^(r,t). 

The aim the paper is to simulate the system with a two-scale semi-lagrangian scheme 

when e ^ 0. Inspired by Frenod, Salvarani and Sonnendriicker [13j, and Frenod, Mouton and Son- 
nendriicker [9], we derive the model (jl.ip by using the two-scale convergence theory developed in 
Allaire [2] and Nguetseng [17J and we obtain a new model which is independent of e. Then, instead 
of discretizing the model with a classical semi-lagrangian scheme, we discretize the new model 
with a semi-lagrangian method in order to obtain an approximation of a function F = F{r^Vr^ r, t), 
where r is the second time scale, which verifies f^{r^Vr^t) ^ 27rF(r, 'U^, Proceeding in such 

a way presents the advantage that there is no longer ^-oscillations in the limit model, so a very 
small time step is no longer required in order to simulate these oscillations. However, since semi- 
lagrangian schemes are based on interpolation on a phase space mesh, we have to pay attention to 
the contributions of the second time scale r in the limit model. For this reason, in the same spirit 
of the work of Lang et a/.[16j, we introduce a r-dependent moving mesh adapted for this two-scale 
numerical method the goal of which is to reduce the number of interpolations during the simulation. 

The paper is organized as follows: in section 2, we recall the procedure leading from the paraxial 
approximation of the complete 3D Vlasov-Maxwell system to the model (jl.ip . and we will recall 
some two-scale convergence results about this model. In section 3, we build the semi-lagrangian 
method on the limit model and we see how to simplify it by considering a particular mesh. Section 
4 is devoted to numerical results obtained with the two-scale semi-lagrangian method: on one 
hand, we compare them to some results obtained from a classical semi-lagrangian method on the 
system in terms of quality of results and CPU time, and on the other hand, we see in this 
section the consequences of the use of the new mesh in the same terms. 

2 The two-scale model 

Firstly, we recall in this section the way to obtain the system from the paraxial approximation 
of the 3D Vlasov-Maxwell equations. Then, we present a theorem about the two-scale convergence 
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of the solutions of (jl.ip which has been proved by Frenod, Salvarani and Sonnendriicker in [13J. 



2.1 Scaling of the paraxial model 

By applying the hypothesis about the considered beam which are mentioned in the introduction, 
the full three-dimensional Vlasov-Maxwell is reduced to 

^.S,/ + v. Vx/ + -(E + H). Vv/ = 0, 

m 

/(x,z = 0,v) =/o(x,v), (2.1) 
-Vx(/> = E, _Ax(/>=^ = -/ fdw, 

where z is the longitudinal position coordinate, x = (x, y) is the transverse position vector, Vz 
is the constant longitudinal speed, v = {vx-,Vy) is the transverse speed vector, / = /(x, z,v) 
is the distribution function of particles whose charge is q and mass is m, = (/)(x, z) is the 
potential linked with the transverse self-consistent electric field E = E(x, z), and H = H(x, z) is 
the transverse external electric field. More details about this derivation can be found in [6j and 
[8]. In this model, we assume that the external electric field is given by 

H(x, z) = -JYo X + Hi (ui y ) X , (2.2) 

where Hq is a positive constant tension. Hi is a /-periodic tension and uoi is a dimensionless real 
constant. 

In the same spirit of Frenod, Salvarani and Sonnendriicker [T3], we build a dimensionless version 
of the system (|2.ip - (|2.2p by introducing the dimensionless variables x',z',v' defined by 



x = Ax', z = L z' ^ w — Vz^'^ (2.3) 

where A is the characteristic transverse radius of the beam and L is the characteristic length of 
the beam. Moreover, we define the dimensionless quantities /q, E', i^Q, by 

/(Ax',Lz',«,v') = ^/'(x',/,v'), E(Ax',Lz') = !^E'(x',^'), 

q"^ XL qL 

/o(Ax',^,v') = ^^/^(x',v'), = flo^o, (2.4) 

q A Li 

(/)(Ax^LzO = -^^\^,z'), Hi{r) = HiH[{r). 

qL 

With these new variables, the system (|2.ip - (|2.2p can be rewritten under the form 



^, _ HpXqL ^, ^ ^ HiXqL ^, f uJiLz' 



Vv'/' = 0, 



'0 

/'(x',^'=0,v') = /^(x',v'), (2-5) 
-Vx'<^' = E', -A^,<P'= [ fdw'. 

Since the beam is supposed to be long and thin, it is natural to take the ratio 

(2.6) 
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Furthermore, as we want to simulate the beam over a large number of periods of the external 
electric field, we also consider the ratio 

i = - <^^^' 

Finally, we suppose that the external electric field is much stronger than the self-consistent electric 
field and that its oscillations in z direction are of the same order as E, so we consider that 

^A^L 1 R[\qL 

5 — - 7 ^ — J- • l^-oj 

Under all these hypothesis the system ()2.5p reduces to 

r(x,v,t = 0) = /o(x,v), 

(2.9) 



where the primed notations for the variables and the initial distribution have been eliminated, and 
where z' has been replaced by t, j' by E' by E^, (\)' by 0^, and by H^. 

We introduce the polar coordinates (r, 6>, v^, ^6>) linked with (x, v) by the relations 

X = r cos ^ , Vr = Vx cos 9 ^ Vy sin ^ , (o ^cW 

y = r sinO ^ vq = Vy cos — Vx sinO . \ • J 

Since the beam is supposed to be axisymmetric, the system does not depend on 0. Furthermore, 
we assume that the initial distribution is concentrated in angular momentum. Then the system 
([2J]) reduces to (frT|) . 



2.2 Two-scale convergence results 

Since the aim of the paper is to develop a two-scale numerical method in order to simulate the 
model (jl.ip , we have to establish that the solution of this model two-scale converges to a function 
F = F(r, Vr^ r, t) in a certain Banach space, and we have to find a system of equations verified by 
F. These results have been proved in [13j and are recalled in the theorem below. 

Theorem 1. We assume that Hq = 1 and that the initial distribution fo of lil.l]) satisfies the 
following properties: 

(i) fo e L1(R2; \r\drdvr) nL^(R2; \r\drdvr) with p > 2, 

(a) fo{r, Vr) > for all (r, Vr) G M^, 

(Hi) / [r'^ ^ v'^) fo{r,Vr)\r\dr dvr < ^oo . 

Then, considering a sequence of solutions (/^,E^) of lil.l]) and extracting a subsequence from it, 
we can say that f^ two-scale converges to F e L^([0,T] x [0, 27r]; L^(R^; \r\drdvr)) and E^ two- 
scale converges toSr e L^([0,T] X [0,27r]; W^'^/^j^R; \r\dr)) . Furth ermore, there exists a function 
G e L^([0,T];L2(R2; \q\dqdur)) such that 

F(r^ r, t) = G( cos(r) r — sin(r) sin(r) r + cos(r) Vr^t) , (2.11) 
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and {G,£r) is solution of 



dtG 



p27r 

/ sm{a)Sr ( cos((t) q + sin((7) cr^ t) da 
Jo 



dqG 



/ sin(cr; 

^0 



27r 



Hi{uJi a) ( cos(cr) q + sin(cr) Ur) da 



dqG 



cos{a)Sr ( cos(<j) q + sin((j) Ur^ cr, t) da 



27T I (iJ ) 

cos(<j) Hi{coi a) ( cos(<j) q + sin(<j) li^) da 



(2.12) 



G((7, ii^, t = 0) = — fo{q, Ur) , 

- dr(r £r{r,T,t)) = [ G( cos(r) r — sin(r) 'U^, sin(r) r + cos(r) 'U^, t) d'Ur , 

?/;/iere /q(cc;i) z5 equal to 1 if uoi e Q, and otherwise. 

Of course, such a result exists for the solution of the system (|2.9p and can be found in [13j. 



3 The two-scale semi-lagrangian method 

In this section, we develop a two-scale semi-lagrangian method in order to approach the solution 
of (jl.ip . in the case where Hq = 1. As it has been suggested in the introduction, the strategy 
is to discretize the model (|2.1ip - (|2.12p in order to obtain a good approximation of F which can be 
used for approaching f^{r,Vr^t) ~ 27rF(r, 'U^, As in the two-scale PIC-method developed by 
Frenod, Salvarani and Sonnendriicker in [13], there is an advantage by proceeding in such a way: 
since there is no longer -^-frequency oscillations in the system (|2.12p . we do not need a very small 
time step for good simulation. In a first time, we recall the basis of a semi-lagrangian method. 
Then we present a motivation for the development of a two-scale semi-lagrangian method through 
the description of a classical backward semi-lagrangian method on the model (jl.ip . Finally, we 
describe the two-scale numerical method itself and we suggest a new mesh in order to simplify it. 

3.1 The semi-lagrangian method 

In this paragraph, we recall the way to discretize the abstract model 

St/(x,t)+U(x,t)-Vx/(x,t)=0 (3.1) 

with a semi-lagrangian method. For this, we have to consider the characteristics of (|3.ip . which 
are the solutions of 

atX(t) =U(X(t),t). (3.2) 
It is an easy game to remark that / is constant along the characteristics, i.e. 



at(/(x(t),t)) =0, 



(3.3) 
(3.4) 



so we can write 

/(X(t;x,8),t) =/(x,8), 

where X(t; x, s) is the solution (|3.2p with the condition X(5) = x. 

This property of / is used in the semi-lagrangian method as follows: assuming that we know 
the value of / at time tn — At on the mesh points (x^)^^o, at, we use the property (|3.4p to say 
that 

/(x„ + At) = /(X(tn - At; ^^,tn^ At),tn - At) , (3.5) 
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so we have to compute the point X(tn — At; x^, + At) first, then compute /(X(tn — At; x^, + 
At),tn — At) by interpolating f{-^tn — At) on the points (xi)^^o, at in order to obtain an 
approximation of f{^i^tn + At). Sonnendriicker et al have suggested in [19] a way to compute a 
second order approximation of X(tn — At;Xi,tn + At): they discretize the equation ()3.2p with a 
finite difference method in order to obtain the fohowing approximation: 

X(tn - At; x^ , tn + At) = x^ - 2 , (3.6) 

where is solution of 

d, = AtU(x, -d„tn). (3.7) 
In many cases, U(-,tn) is only known at points x^. Then we have to replace (|3.7p by 

di = At nU(xi - di, tn) , (3.8) 

where 11 is an interpolation operator on points (xi)^. Assuming that the polynomial function 
X ^ nU(x, tn) is regular enough, we write the following expansion of ()3.8p : 

d, = AtU(x,,tn) - At Vx(nU)(x,,tn) d, + 0{At^) , (3.9) 

because, in the expansion of IIU in x, we get a (9(|dip) term, which is a (9(At^). Then, we obtain 
a second order accurate approximation of d^ given by 

d, = At (Id + At Vx(nU)(x,,tn))"' X U(x„tn) . (3.10) 
Considering this approach, the semi-lagrangian method writes: 

1. knowing / at time tn — At and U at time t^, we compute (di)i = o, at by using the relation 
(|3.1()p for each i, 

2. we compute / at time tn + At as follows: 

/(x„ tn + At) = n/(x, - 2 d,, tn - At) . (3.11) 

3.2 Implementation of the non-homogenized model 

In this paragraph, we describe a classical semi-lagrangian method on the system (|l.ip . Since the 
electric fields and do not depend on Vrj we can do a time-splitting on the first equation of 
the model, i.e. solving separately at each time step the equations 

dtr + ^drr=0, (3.12) 

and 

dtr^{E^r^El)d,J' = 0, (3.13) 

with a second order in time numerical scheme instead of solving the complete equation with the 
same scheme. As a consequence, we only do ID interpolations instead of 2D interpolations. Then, 
denoting with the aproximation of /^(•, •, t^) and with the approximation of E^{-^ tn) on the 
uniform mesh {ri^Vrj)ij where Ar is the size of one cell in r direction and A^'^ is the size of one 
cell in Vr direction, an iteration of the method is organized as follows: 

1. Knowing and E^, we do a backward advection of Ai in Vr direction and we define /| by 

mn,vrj)=n,j:,{n,vr, - ^{E^Jn) + El{n,tn))) , (3.14) 

where 11^^ is a ID cubic spline interpolation operator on the points {vrj)j- 
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2. We do an advection of At in r direction and we define by 

where 11^ is a ID cubic spHne interpolation operator on the points (r^)^. 

3. We compute E^_^^ by discretizing the formula 

1 



(3.15) 



'i Jo 



s fl^{s,Vr)dvrds if Vi ^ 0, 
otherwise, 



(3.16) 



with the trapezoidal rule on the points {ri,Vrj)ij. 
4. We compute by doing a last advection of ^ in Vr direction: 



(3.17) 



If we use such a method, we have to guarantee the accuracy of the scheme, especially we 
consider the ^-frequency oscillations of the external electric field. A solution is to assume that the 
time step At satisfies 



Ar < 



At 



At 



(3.18) 



Vrj 



AVr < Vrj - —{E^{ri) ^E^{ri,tn)) < Vrj^AVr, 



for all i^jjTi. Then, in order to obtain good results with this method. At has to be of the same 
order of e which penalizes the method in terms of CPU time cost when we consider a very small e. 



3.3 Implementation of the two-scale model 

In this paragraph, we adapt the semi-lagrangian method we have described in the paragraph 3.1 
to the model ()2.12p . In this case, the characteristics of the system are the solutions of 



dtQ{t) = {Si){Q{t),Ur{t),t) 

dtUrit) = {S2){Q{t),Ur{t),t) 



(3.19) 



where (Si) and {£2) are defined by 



{Si){q,Ur,t) 



{S2){q,Ur,t) 



sin(cr) £r{cos{a) q -\- sm{a) Ur, cr^t) 

_^ Iq{^i) / ^.Qgj^^j ^ ^ sin(cr) Ur) 

ZTT 

f 

/ cos(cr) Sr{cos{a) q -\- sm{a) Urj cr^t) 
I ( ) 

+ Hi{(jOi (j) ( cos((t) q + sin((j) Ur) 



da . 



da . 



(3.20) 



As in the paragraph 3.1, we remark that the solution G of (|2.12p is constant along the character- 
istics, so we can write: 



G(g, Ur^ ^n+l) — G{Q{tn-l]q-> ^n+l), Ur{tn-l]Ur^ ^n+l), ^n- 



(3.21) 



where tn = n At, and where (Q(tn-i; ^, ^n+i), ^r(^n-i; '^r, ^n+i)) is the solution of (|3.19p with the 
condition (Q(tn+i), /7r(tn+i)) = {q,Ur). 
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Firstly, we define the mesh in q and Ur directions by considering the points Qi = i Aq and 
Urj = j Aur {i = —Pq, . . . ^ Pq, j = —Pur^^ • • • , Pu^) • Wc also considcr the uniform mesh Tm = m At 
on [0, 27r] (m = 0, . . . , Pr). Finahy, we fix a time step At for the entire simulation. Then, denoting 

with the approximation of G at time an iteration of the semi-lagrangian method is organized 
as follows: 

1. Assuming that we know the value of and G^~^ on the mesh {qi,Urj)ij, we compute 
with the formula 



Qi Jo JR 

sin(r^)s + cos(r^)'U^) ds dvr if i ^ 0, 



(3.22) 







otherwise. 



Since G^ is only known at points {qi^Urj), we have to interpolate G^. Assuming that the 
support of is included in [-{Pq + l)Aq, {Pq + l)Aq] x [-(Pn, + l)Aur, {Pu^ + l)Aur\ 
with Pg, Pu^ G N large enough, we use the trapezoidal rule in order to approach the integral 
above. We obtain 



Aq Aur 



^ n2G''(cos(r,n)<7i - ^in{Tm)Urj,^in{Tm)qi + cos(r^)? 



3 = -Pur. 

i-1 



(3.23) 



- ^ kU2G'^{cos{rm)qk - sm{rm)urj,sm{rm)qk + cos(r^)i 



where 112 is a cubic spline interpolation operator on the points (qi^Urj). 

2. We compute (ff) and {£2) at points {qi^Urj): since 8^ is only known at points (qi^Tm), 
we have to interpolate By using the trapezoidal rule for approximating the integrals in 
(|3?2Q|) . we obtain 



{S^){qi,Urj) 



{£^){qi,Urj) 



-At ^ sin(rm) IIi^;? ( cos(r^) + sin(r,n) Urj , T77 



m = 



27r 



i^i(^i r,n) ( cos(r^) qi + sin(r,n) ^ 



Ar ^ cos(r^) Ilif;:' ( cos(r^) + sin(r^) 

^q(^i) rr / \ / / \ , • / \ 

■ Hi[u;i Tm) [ cos(r^) + sm(r^) ^ 



m = 



where Hi is a cubic spline interpolation operator on the points {qi)i. 
3. We compute the the shifts {dj.^d'^ ,) by using the following formula: 



4i 



where the matrix Aij is defined by 



Aij = Id + At 



dq {U2 {£l )){qi, Urj ) dur. (f ^ ) ) ((7i , ^^rj ) 

4. We compute G^+-^ by interpolating G^~^ at the points {q^ — 2dlj^Urj — .,): 
G^^\q,,Urj) = U2G^-\q, -2dl^,Urj - 24^.) . 



(3.24) 



(3.25) 



(3.26) 



(3.27) 
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5. We save the approximation of at time t^+i given by 

r (r, Vr.tn^i) - 2 TT EsG^+i ( COS (^) r-sin (^) sin (^) r+cos (^) . (3.28) 

In order to initialize this two time step advance, we have to compute from which is given 
as an initial data: for this purpose, we perform a complete iteration such as decribed above where 
At is replaced by ^ and where we assume that G^/^ = G^ . 

3.4 The two-scale mesh 

In most of test cases, we assume that the support of the initial distribution /o is compact and is 
included in some n = [-R, R] X [-VR, Vr] C for R>0,vr>0 large enough. Then, if we follow 
the algorithm presented in the previous paragraph, the first thing we have to compute is by 
approximating the integral in (|3.22p : for numerical reasons, we have to reduce the integral on R to 
an integral on a compact interval. Furthermore, since we can only say that the support of (r, Vr) ^ 
/o ( cos(r) r— sin(r) Vr^ sin(r) r+cos(r) Vr) is included in = [—R—vr^ R-\-vr] x [—R—vr^ R-\-vr] as 
it is illustrated in Figure 1, the integral on R in (|3.22p is reduced to an integral on [—R — vr^ R-\-vr]. 
As a consequence, we have to do all the simulation on ft' instead of ft in order to avoid losing 
some data and we have to increase the number of mesh points in order to keep good interpolation 
results, even if the distribution function will be equal to at many new points. 



Vr 



VR^R 



Vr 



-R-vr 



-R 



R 



R^VR 



-Vr 



-Vr - R 



Figure 1: Support of (r, Vr) ^ /o(^, Vr) and (r, Vr) ^ /o ( cos(r) r — sin(r) Vr, sin(r) r + cos(r) Vr) 
with r = ^. 



In this paragraph, we present a different approach in order to avoid this extension of the 
simulation domain and its mesh. Before explaining the main idea of this new method, we consider 
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the following meshes on Q and [0,27r]: 



M{Q) = {(r^, Vrj) = (i Ar, j Avr) : i = -Pr, ...,PrJ = -P^,, . • • , } , 

M([0,27r]) = ^Tm = mAT : m = 0,...,P^}, (3.29) 
M([-i?,i^]) = {r^ =2Ar : i = -P^, . . . , P^} , 



where Ar = pqij , = p^^^ and Ar = • Considering the function 7 defined by 



7 : M2 X 27r] 
(r, '?;^,r) 



(cos(r) r — sin(r) t'^-, sin(r) r + cos(r) Vr) 



(3.30) 



we define ^}{t) and M(r^(T)) by 



n(r) =7(n X {r}) CM^ 
M(l7(r)) = {7(n,^^,-,T) : i = -P^, . . . , P^ , j = -P^^, . . . , P^ J 



(3.31) 



^ ^_ 



* i ♦ 



4. ^ ^ HH 



4. ^_ 



- + + - 



Figure 2: Mesh M(0) = M(0(0)) and support of {r,Vr) ^ fo(r,Vr). 
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Figure 3: Mesh M(Q{^)) and support of (r, '^r) ^ /o o 7(^, '^r, f )• 



Figure 4: Mesh M(Q{^)) and support of (r, Vr) ^ /o o 7(^, ^r, 

The main idea of our new method is to compute (r^Vr) ^ F{r^Vr^Tm^tn) at points {vi^Vrj) G 
M(ll) and to compute r fr(^, ^m^tn) at points G M([— i?, i?]), whereas the function (g', Ur) ^ 
G{q^Ur,tn) is computed at points ^{ri^Vrj^Tm) G M(l](rm)) for every G M([0,27r]). This 
approach is similar as the time-dependent moving grid described by Lang et al. in [16], even if, in 
our case, the mesh M(r^(r)) only depends on M{VL) and r, and is completely defined before the 
beginning of the simulation. 
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As a first consequence, considering that the support of {r^Vr) ^ /o(r, v^) is included in Q is 
equivalent to considering that the support of (r^Vr) ^ /o(7(^, '^r, 't)) is included in ^{r) for any 
r G [0, 27r] such as illustrated in Figures 2,3, and 4, where the support of a Kapchinsky-Vladimirsky 
distribution is represented (see [15], [18] or [8] for more details about this distribution). Then we 
do not have to extend Q in order to avoid losing some data. Furthermore, the equation (|2.1ip 
reads 

F{r,Vr,T,t)=G{-f{r,Vr,T),t) . (3.32) 

Considering this approach, we fix a time step At for all the simulation. Then, an iteration of 
this new semi-lagrangian method is organized as follows: 

1. Assuming that G^~^ and are known on the mesh M(Q{rm)) for all G M([0,27r]), 
we compute at points {vi^Tm) G M(^[—R,R]) x M([0,27r]). With the new notations, the 
equation ()3.22p simplifies itself to 



— / / sG''{-f{s,Vr,rm)) dvrds ifiy^' 



(3.33) 



otherwise. 



Assuming that the support of (7(-, •, r^n)) is included in l^(rm) for each (which is 
equivalent to assume that the support of G"^ (7(-, •, 0)) is included in ft), we use the trapezoidal 
rule to approximate the integral above: 



Ar Avr 



i-l 



G" {iin,vr,, Trn)) + fcG"(7(rfc,«,„T™)) . (3.34) 



k = l 



We remark here that, contrary to the computation done in (|3.23p . we do not have to inter- 
polate G^. 

2. We compute and {E2) at points ^{ri^Vrj^Tm) ^ M{p.{Tm))'' 



27r 



Sin(cr) f^(cOs(cr - Tm) Vi + sin(cr - T^) Vrj.CT 



27r 



Hi{u;i a) ( cos(cr - r^) + sin(cr - r^) v. 
Jo 

' Hi{u;i a) (cos(cr - rm)ri + sin(cr - r^) Vr 



da . 



da . 



(3.35) 

We approximate the integrals with the trapezoidal rule. However, since £^ is only known 
at points of M{[—R^R]) x M([0,27r]), we have to interpolate it: for that we choose a cubic 
spline interpolation operator on the mesh M(^[—R,R]) and we denote it with IIi. Then we 
have the following approximations: 

Pr 

{^l){l{U,Vrj,rm)) ^-Ar^ Sin(r/e) XIi^;? ( COs(r/e -Tm)ri + Sin(r/e -Tm)Vrj,Tk) 



k = 



27r 



^1(^1 Tk) (cos(rfe - r^) + sin(r/e -rm)n 



{£^){j{ri,Vrj,rm)) ^Ar^ cos(r/e) UiS^!' {cos{Tk -rm)ri + sin(rfe - r^) ^rj 



k = 



/q(cji) 
27r 



^1(^1 Tk) ( cos(rfe - r^) n + sin(r/e - Tm) Vrj 



(3.36) 



12 



3. We compute the shifts d{ri,Vrj,Tm) verifying 

For that, we consider the cubic sphne interpolation operator 112^ on the mesh M(^}{rm)) 
and, inspired by (|3.1Qp . we consider the fohowing approximation of d(ri, v^^, r^n): 

-1 f {^l){7{ri,Vrj,rm)) 



d(n,.,,,.„, = A*A-„^;-;;;;;;;.;;~j J, (3.38) 

where the matrix Aij^rn is defined by 
Ai^j^m = Id + AH . (3.39) 

4. We compute G^+^ on the meshes M{Q{rm)) for ah G M([0,27r]): 

G^+i (7(r„ Vrj.Tm)) = U^G^-' (7(r„ r^) - 2 d(r„ , r^)) (3.40) 

5. Assuming that there exists a fixed integer K eN such that 

At = eArK, (3.41) 
we save the approximation of on M{^}) given by 

r{n,Vrj,tn+i) ^ 2^G"+^(7(r„^.„r(,+i)K)) • (3.42) 

Contrary to the scheme described in the previous paragraph, we do not have to interpolate 
G"^+^ to obtain an approximation of at time t^+i- 

Concerning the initialization of this two time step advance, we compute = ^ fo on the meshes 
M(Q{rm)) for all G M([0,27r]) and, assuming that G^/^ = G^, we compute G^ by performing 
a complete iteration such as described above with At replaced by 

We have built another two-scale semi-lagrangian method based on a mesh depending on the 
variable r. Compared to the semi-lagrangian method described in the paragraph 3.3, this technique 
allows us to avoid some interpolations not only when we compute the two-scale limit of the electric 
field S but also when we build the approximation of defined by (|3.42p . As a consequence, we 
reduces the global numerical diffusion within the simulation. However, we have to compute G and 
{£) on each mesh M(Q{rm))j which can be expensive in CPU time. 

4 Numerical results 

In this section, we present some numerical results obtained with the two-scale semi-lagrangian 
methods described in the paragraphs 3.3 and 3.4. Following the approach of |T3] for validating our 
new methods, we study in a first time some linear cases where we can find an analytic expression 
of the solution G of the system (|2.12p . Then we test the methods on non-linear cases. 

4.1 Linear cases 

In order to validate the two-scale semi-lagrangian methods described in the paragraphs 3.3 and 
3.4, we simulate some linear cases, i.e. with a self-consistent electric field set to 0. This assumption 
allows us to compute analytically the solution G of the system (|2.12p under an adequate choice of 
LOi and Hi. 
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t = 1.1088 t = 6.468 

Figure 5: Simulations of type (I) (first row), (II) (second row), (III) (third row), and (IV) (fourth 
row) for a semi-gaussian beam wihout self-consistent electric field and coi = 4>/2, i^i(r) = cos(r). 

In a first example, we suppose that coi ^ Q, so we obtain that G is stationary in t, and is equal 
to ^ /o- Then, the two-scale simulation reduces itself to the computation of 

(r, t) I — > /o ^ cos (-) r — sin (-) sin (-) r + cos (-) , (4.1) 

for any function Hi . It is much simpler than simulating the model (jl.ip with the semi-lagrangian 
method described in the paragraph 3.2. In Figure 5, we observe such a case, with uji = 4>/2, 
Hi{t) = cos(r), and /o given by 

2 

fo{r,Vr,t) = ^ exp ( - X[-r^,r^](0 , (4-2) 

V 27r vth ^ ^^th"^ 

with X[-r^,r^](^) = 1 if 1^1 < ciud othcrwisc. This corresponds to a semi-gaussian beam in 
particle accelerator physics. In these Figures, we suppose that Tm = 0.75, Vth — 0.1, no = 4 and 
e = 10~^. Furthermore, the simulations (I), (II), (III) and (IV) correspond to 
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• simulation (I): we solve the system (jl.ip with a classical semi-lagrangian method, with Pr = 
Py^ = 64 and R = vr = 3, 

• simulation (II): we solve the system with a classical semi-lagrangian method, with 
Pr = Py^ = 128 and R = vr = 3, 

• simulation (III): we solve the system (|2.12p with a two-scale semi-lagrangian method on a 
two-scale mesh, with Pr = Py^ = 64, P^- = 16 and R = vr = 3, 

• simulation (IV): we solve the system ()2.12p with a two-scale semi-lagrangian method on a 
uniform mesh, with Pq = P^^ = 128, Pr = 16, and R = vr = 3, and we compute the 
approximation (|3.28p on a uniform 129 x 129 grid in (r, v^) on [—R^R] x [—vr^vr]. 

Since we can compute an analytic solution G of (|2.12p . we can compare the approximations of 
the solution of (jl.ip to the function defined by ()4.ip . These comparisons are summarized in 
Figure 6: in this figure, we present the norm of the difference between the function (14.11) and 




time 



Figure 6: Evolution of the norm of the difference between the function ()4.ip and the approxi- 
mation fh computed with the simulations (I), (II), (III), and (IV). 

In a second example, we suppose that uj G N>2. Then, if we assume that i^i(r) = cos^(r), the 
model (|2.12p reduces itself to 

dtG-^dqG^^du^G = 0. (4.3) 
Knowing the initial data /o, we can write 

G((7, Ur,t) = ^ fo(^ cos(^) q + sin(^) Ur, - sin(^) q + cos(^) Ur^ , (4.4) 

so the two-scale simulation reduces itself to the computation of 

{r,Vr,t) I — >/o(^cos(^ - ^)r + sin(^ - ^) v^,-sin(^ - ^)r + cos(^ - ^) ^r) • (4.5) 

In Figure 7, we observe such a case with coi =2, /o given by (|4.2p with Tm = 0.75, Vth = 0.1 and 
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t = 0.2957 t = 5.9875 



Figure 7: Simulations of type (I) (first row), (II) (second row), (III) (third row), and (IV) (fourth 
row) for a semi-gaussian beam wihout self-consistent electric field and cji = 2, i^i(r) = cos^(r). 

As we can see in Figures 5 and 7, the classical semi-lagrangian method needs a very refined 
mesh in r and Vr directions in order to produce good results. Otherwise, it produces results which 
do not correspond to physics. We can explain it by the small time step induced by the condi- 
tion (|3.18p in order to guarantee the robustness of the method, and then a very high number of 
interpolations which introduce numerical diffusion. As an example, we have taken a time step 
AtNH ~ 1.5 X 10~^ for the simulation of type (I) of the second case. 

On the other hand, both two-scale numerical methods we have described in paragraphs 3.3 
and 3.4 produce good results, even if we consider a non-refined mesh in r and Vr directions: this 
phenomenon can be explained by the independence of the time step in e. Consequently, we are 
allowed to take a time step Atn much larger than we can in the classical semi-lagrangian context, 
and then, we can significantly reduce the number of interpolations and the numerical diffusion 
which is induced. For example, for simulations of type (III) and (IV) on the second test case, we 
have defined Atn with the formula ()3.4ip where K = 2, giving us Atn ^ 7.4 x 10~^. 

Finally, by observing Figure 6, we remark that the use of the two-scale mesh in the first case 
reduces significantly the error between the function defined by (|4.ip and the approximation of 
given by the discretization of the two-scale model ()2.12p : indeed, the norm of this error is 
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nearby 10~^ when the approximation of is given by the simulation (III) whereas it oscihates 
between and 0.75 when the approximation of is given by the simulation (IV). Furthermore, 
the oscillations of this error for simulations (I) , (II) , and (IV) are due to the fact that we discretize 
a non-smooth semi-gaussian distribution the support of which rotates in the phase space under the 
action of the external electric field and that we do it on a uniform phase space mesh. 

As a first conclusion, we can say that, on linear cases, the two-scale semi-lagrangian methods 
we have proposed give much better results than the classical semi-lagrangian method on the same 
mesh in r and v^, which is promising for non-linear cases. Furthermore, as it was announced in the 
paragraph 3.4, the use of the two-scale mesh reduces significantly the numerical diffusion linked 
with the global number of interpolations within the simulation. 

4.2 Non-linear cases 




t = 1.4784 t = 3.234 t = 5.544 



Figure 8: Simulations of type (IF) (first row), (III) (second row), and (IV) (third row) for a 
semi-gaussian beam with i^i(r) = cos(r) and coi = 4>/2. 

In this paragraph, we do not assume that the self-consistent electric field vanishes. Then, in 
most of cases, we are not able to find an analytic expression of the solution G of (|2.12p . So, in 
order to validate our two-scale methods, we have to compare their results on such a case to the 
results produced by a classical semi-lagrangian method on (|l.ip with the same initial data, and 
pay attention to the development of thin structures within the beam. For that, we consider the 
simulations (I), (IF), (III) and (IV), where the simulation type (IF) corresponds to a classical 
semi-lagrangian method on the system (|l.ip with Pr = Py^ = 256 and r = Vr = 3. Since the 
simulation (I) already gives bad results on linear cases, it is not really useful to present its results 
in terms of quality on non-linear cases. However, we have to pay attention to its CPU time cost 
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in order to compare it to the other simulation types one^. 

In a first case, we suppose that i^i(r) = cos(r), c^i = 4>/2, e = 10~^, and /o is given by (|4.2p 
with Tm = 0.75 and Vth = 0.1. In order to guarantee the robustness of the two-scale schemes, we 
suppose that the time step Atn is computed by using the the formula ()3.4ip with K = 5, giving 
us AtH ^ 0.0185, whereas the time step AtNH for the simulations (I) and (IF) is of the form 
with N large enough in order to verify the condition (|3.18p . In Figure 8, we observe the results 
obtained with simulations of type (IF), (III), and (IV), and we can find some results in terms of 
CPU time in the Table 1. 

In a second case, we suppose that Hi{r) = cos^(r), cji = 2, e = 10~^, and /o is given by (|4.2p 
with Tm = 0.75 and Vth = 0.1. In order to guarantee the robustness of the two-scale schemes, 
we suppose that Atn is defined by (|3.4ip with K = 2, which gives us Atn ^ 7.392 x 10~^, and 
AtNH = with N large enough. In Figure 9, we observe the results obtained with simulations 
of type (IF), (III), and (IV) and CPU times costs for each simulation can be found in the Table 1. 




t = 1.1458 t = 3.6221 t = 5.8027 



Figure 9: Simulations of type (IF) (first row), (III) (second row), and (IV) (third row) for a 
semi-gaussian beam with cji = 2, Hi{r) = cos^(r). 

As we can see in Figures 8 and 9, the three simulations (IF), (III) and (IV) produce results 
of the same quality, even if the mesh in r and Vr directions used in the classical semi-lagrangian 
method is much more refined than the one used for the two-scale simulations, which enlarge to 
non-linear cases the conclusion we have established at the end of the previous paragraph. 



^All the simulations have been conducted on a Sun Fire X4600 server with an AMD Opteron 8220 processor 
(Dual Core 3000 Mhz) under SunOS 5.10 system. 
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Case 


Simulation (I) 


Simulation (IF) 


Simulation (III) 


Simulation (IV) 


CPU time 


TV 


CPU time 


TV 


CPU time 


CPU time 


iJi — 4>/2, Hi — cos 


35m 


122 


35h 6m 50s 


480 


Ih 43m 39s 


55m 3s 


LJi = 2, i^i = cos'^ 


37m 32s 


49 


38h 7m 6s 


192 


5h 45m 25s 


2h 37m 25s 



Table 1: CPU time costs: the final time is T = 6.93 for the case where uji = 4V^ and Hi = cos, 
and T = 6.9854 for the case where uoi = 2 and Hi = cos^. 



Furthermore if we observe the CPU times costs of each simulations (see Table 1), we remark 
that both two-scale numerical methods are much slower than the classical semi-lagrangian method 
when they are ran on a same mesh in r and Vr directions. It is not surprising because of the high 
number of fixed point problems we have to solve during the two-scale simulations. But on the other 
hand, if we compare the CPU time costs of the simulations (IF), (III) and (IV) which give the 
same quality of results, we remark that both two-scale numerical methods we have described are 
much faster than a precise classical semi-lagrangian method. One more time, this phenomenon can 
be explained by the condition (|3.18p imposed to the time step AInh within the classical method: 
if we refine the mesh of the phase space, we diminish the time step, and then increase the number 
of iterations in time we need to reach the the final time of the simulation. As a conclusion, we can 
say that if we want high quality results, it is preferable in terms of CPU time cost to run one of 
the two-scale methods we have developed instead of the classical method. 

In a last case, we suppose that Hi{r) = cos^(r), cji = 1, e = 10~^, and /o is given by (|4.2p 
where no = 4, r^n = 1-85 and Vth = 0.1. As in the previous tests, we suppose that the time step 
Atn for the two-scale methods is given by (|3.4ip where K = 1, and that the time step AInh for 
the classical semi-lagrangian method is of the form with N large enough. The goal of this 
numerical experiment is to observe the same structures as Frenod, Salvarani and Sonnendriicker 
have observed in [13j. Since the structures we want to observe are quite thin, we consider the 
simulation types (IF'), (IIF) and (IV) corresponding to 

• simulation (IF'): we solve the system (jl.ip with a classical semi-lagrangian method, with 
Pr = Pvr. = 256 and R = vr = 3, 

• simulation (IIF): we solve the system (|2.12p with a two-scale semi-lagrangian method on a 
two-scale mesh, with Pr = Py^ = 128, Pr = 20 and R = vr = 3, 

• simulation (IV): we solve the system (|2.12p with a two-scale semi-lagrangian method on a 
uniform mesh, with Pq = P^^ = 256, Pr = 20, and R = vr = 3, and we compute the 
approximation (|3.28p on a uniform 257 x 257 grid in (r, v^) on [— P, P] x [—vr^vr]. 

In the Figure 10, we can observe some results obtained with the simulations (IF'), (IIF) and 
(IV). One more time, we can remark that the two-scale methods do not need a mesh as refined 
as the classical method's one for producing high quality results. Furthermore, we can observe 
in the last column of the Figure 10 that the beam simulated with the classical semi-lagrangian 
method becomes unfocused in long time, even if we consider a highly refined mesh. This remark 
is confirmed by the Figure 11 where we notice the thin structures are better determined with a 
two-scale method. Moreover, the result given by the classical method is slightly out of phase and 
more diffusive. One more time, the main reason of this problem is the condition (|3.18p imposed on 
the time step in the classical semi-lagragian simulation: if we consider a 513 x 513 grid in (r, v^), 
this condition induces a so small value of AInh (nearby 3.84 x 10~^) that the numerical diffusion 
introduced by the interpolations makes the beam unfocused in a long time simulation. On the 
other hand, the two-scale results match well with the expected long time behavior as described in 

m 
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t = 1.3464 t = 4.3388 t = 5.1462 



Figure 10: Simulations of type (11") (first row), (HI') (second row) and (IV) (third row) for a 




(II") (iir) 

Figure 11: Simulations of type (11") and (HI') for a semi-gaussian beam at time t = 5.984 with 
cji = 1, Hi{t) = cos^(r) and = 1.85. 

5 Conclusion and perspectives 

We have built a two-scale semi-lagrangian method and proposed a new mesh in order to simplify 
the computation of the electric field and, which leads to another two-scale semi-lagrangian method. 
These methods have been tested on non-smooth initial conditions, especially semi-gaussian beam 
initial conditions: on linear cases, we have concluded that these two-scale methods are very efficient. 
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even if we consider a coarse mesh in r and Vr directions, contrary to a classical semi-lagrangian 
scheme on the model (jl.ip which needs much more points to produce good results. On non- 
linear cases, we have reached the same conclusion for short time simulation, not only in terms 
of quality of results, by also in terms of CPU time costs. These results are very promising for 
extensions to higher dimensional problems such as the two-scale limit models obtained by Frenod 
and Sonnendriicker in [iO\ \TT\ \12\ , the finite Larmor radius approximation obtained in [4J, or other 
charged particle beam problems which cannot allow any time splitting. Furthermore, for long time 
simulation, both two-scale numerical methods we have developed give very good results, contrary 
to the classical semi-lagrangian method which is so penalized by its numerical diffusion for long 
time simulations that it produces results which do not not correspond to the expected behavior. 
These results are also promising since they consolidate the conclusions of Frenod, Salvarani and 
Sonnendriicker in [13], which are that a two-scale numerical method can be successfully used in a 
context of non-smooth initial data. 

Finally, we also have remarked that, even if they are faster than the classical semi-lagrangian 
method for obtaining the same quality of results, both two-scale numerical methods need a very 
high CPU time cost. Since this time is essentially spent for computing the fixed point problems 
within the methods, it can be interesting to find a way to improve this part of the methods. 

Acknowledgements The author wishes to thank E. Frenod and E. Sonnendriicker for the dis- 
cussions on the topic. 
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